{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation\n",
    "\n",
    "## The dice problem\n",
    "\n",
    "Suppose I have a box of dice that contains a 4-sided die, a 6-sided die,\n",
    "an 8-sided die, a 12-sided die, and a 20-sided die. If you have ever\n",
    "played *Dungeons & Dragons*, you know what I am talking\n",
    "about.\n",
    "\n",
    "Suppose I select a die from the box at random, roll it, and get a 6.\n",
    "What is the probability that I rolled each die?\n",
    "\n",
    "Let me suggest a three-step strategy for approaching a problem like\n",
    "this.\n",
    "\n",
    "1.  Choose a representation for the hypotheses.\n",
    "\n",
    "2.  Choose a representation for the data.\n",
    "\n",
    "3.  Write the likelihood function.\n",
    "\n",
    "In previous examples I used strings to represent hypotheses and data,\n",
    "but for the die problem I’ll use numbers. Specifically, I’ll use the\n",
    "integers 4, 6, 8, 12, and 20 to represent hypotheses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from dice import Dice\n",
    "\n",
    "suite = Dice([4, 6, 8, 12, 20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And integers from 1 to 20 for the data. These representations make it\n",
    "easy to write the likelihood function:\n",
    "\n",
    "```python\n",
    "class Dice(Suite):\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if hypo < data:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1.0/hypo\n",
    "```\n",
    "\n",
    "Here’s how `Likelihood` works. If `hypo<data`, that means the roll is\n",
    "greater than the number of sides on the die. That can’t happen, so the\n",
    "likelihood is 0.\n",
    "\n",
    "Otherwise the question is, “Given that there are `hypo`\n",
    "sides, what is the chance of rolling `data`?” The answer is\n",
    "`1/hypo`, regardless of `data`.\n",
    "\n",
    "Here is the statement that does the update (if I roll a 6):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.08500000000000002"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "suite.Update(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here is the posterior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 0.0\n",
      "6 0.3921568627450979\n",
      "8 0.2941176470588235\n",
      "12 0.19607843137254896\n",
      "20 0.11764705882352941\n"
     ]
    }
   ],
   "source": [
    "suite.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we roll a 6, the probability for the 4-sided die is 0. The most\n",
    "likely alternative is the 6-sided die, but there is still almost a 12%\n",
    "chance for the 20-sided die.\n",
    "\n",
    "What if we roll a few more times and get 6, 8, 7, 7, 5, and 4?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for roll in [6, 8, 7, 7, 5, 4]:\n",
    "    suite.Update(roll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this data the 6-sided die is eliminated, and the 8-sided die seems\n",
    "quite likely. Here are the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 0.0\n",
      "6 0.0\n",
      "8 0.9432484536722127\n",
      "12 0.055206128061290875\n",
      "20 0.001545418266496554\n"
     ]
    }
   ],
   "source": [
    "suite.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the probability is 94% that we are rolling the 8-sided die, and less\n",
    "than 1% for the 20-sided die.\n",
    "\n",
    "The dice problem is based on an example I saw in Sanjoy Mahajan’s class\n",
    "on Bayesian inference. You can download the code in this section from\n",
    "<http://thinkbayes.com/dice.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The locomotive problem\n",
    "\n",
    "I found the locomotive problem in Frederick Mosteller’s, *Fifty\n",
    "Challenging Problems in Probability with Solutions* (Dover,\n",
    "1987):\n",
    "\n",
    "> “A railroad numbers its locomotives in order 1..N. One day you see a\n",
    "> locomotive with the number 60. Estimate how many locomotives the\n",
    "> railroad has.”\n",
    "\n",
    "Based on this observation, we know the railroad has 60 or more\n",
    "locomotives. But how many more? To apply Bayesian reasoning, we can\n",
    "break this problem into two steps:\n",
    "\n",
    "1.  What did we know about $N$ before we saw the data?\n",
    "\n",
    "2.  For any given value of $N$, what is the likelihood of seeing the\n",
    "    data (a locomotive with number 60)?\n",
    "\n",
    "The answer to the first question is the prior. The answer to the second\n",
    "is the likelihood.\n",
    "\n",
    "![Posterior distribution for the locomotive problem, based on a uniform\n",
    "prior.](figs/train1.pdf)\n",
    "\n",
    "[fig.train1]\n",
    "\n",
    "We don’t have much basis to choose a prior, but we can start with\n",
    "something simple and then consider alternatives. Let’s assume that $N$\n",
    "is equally likely to be any value from 1 to 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hypos = range(1, 1001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now all we need is a likelihood function. In a hypothetical fleet of $N$\n",
    "locomotives, what is the probability that we would see number 60? If we\n",
    "assume that there is only one train-operating company (or only one we\n",
    "care about) and that we are equally likely to see any of its\n",
    "locomotives, then the chance of seeing any particular locomotive is\n",
    "$1/N$.\n",
    "\n",
    "Here’s the likelihood function:\n",
    "\n",
    "```python\n",
    "class Train(Suite):\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if hypo < data:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1.0/hypo\n",
    "```\n",
    "\n",
    "This might look familiar; the likelihood functions for the locomotive\n",
    "problem and the dice problem are identical.\n",
    "\n",
    "Here’s the update:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxQAAAEtCAYAAACRaSbdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XtcVHX+P/DXzMDADPeLDFdBRUXBEhlMRA0vQK5t62XV\nci3NtWKzC9jPymptK+m7tq7Xr4K7mmm6ZV9zrXbVQE1DjRQMEgGvgHIZLnK/CQzz+wM5MgJyZ2B4\nPR8PHp458zln3odOyovP53M+Io1GowEREREREVEniHVdABERERER9V8MFERERERE1GkMFERERERE\n1GkMFERERERE1GkMFERERERE1GkMFERERERE1GkMFERERERE1Gl9IlBs374dQ4cOhUwmg1KpxJkz\nZx7aPikpCQEBAZDL5XBxccFHH33UrM3p06ehVCohk8ng7u6OHTt2NGtTVlaG1157DU5OTjA2NsaI\nESNw8ODBbrsuIiIiIiJ9Z6DrAg4cOIDQ0FBERkbC398f27Ztw8yZM5GSkgJnZ+dm7cvKyhAYGIiA\ngADEx8cjJSUFS5cuhampKcLCwgAA6enpmDVrFpYvX479+/cjJiYGL7/8Muzs7DBnzhwAQF1dHWbM\nmAFbW1scPHgQTk5OyMzMhJGRUa9ePxERERFRfybS9UrZEyZMwNixYxEZGSnsGzFiBObPn4/w8PBm\n7SMiIrB69Wrk5eVBKpUCAMLDwxEZGYnbt28DAN566y0cPnwYV65cEY574YUXkJycjLNnzwIA/vGP\nf+CTTz5BamoqDAx0nquIiIiIiPolnQ55qq2tRXx8PAIDA7X2BwUF4dy5cy0eExsbi8mTJwthAgCC\ng4ORnZ2NjIwMoU1QUJDWccHBwYiLi4NarQYAfPPNN/D398crr7wCBwcHeHp64oMPPkBdXV13XiIR\nERERkV7TaaAoKCiAWq2GQqHQ2q9QKKBSqVo8RqVStdheo9EIx7TWpq6uDgUFBQCAmzdv4quvvkJd\nXR2OHDmCtWvXIjIyEu+88053XR4RERERkd4bsGN96uvroVAo8M9//hMikQje3t4oKCjAypUr8ckn\nn2i1LSkp0VGVRERERERdZ2Fh0WPn1mmgsLW1hUQiQW5urtb+3Nxc2Nvbt3iMvb19i+1FIpFwTGtt\nDAwMYGtrCwBwcHCAVCqFSCQS2owaNQqVlZW4c+cObGxsunx9RERERET6TqdDngwNDeHj44Po6Git\n/dHR0fD392/xGD8/P8TExKCmpkbYFxUVBUdHR7i6ugptHjxnVFQUlEolJBIJAMDf3x/Xr1/XanPl\nyhXI5XKGCSIiIiKidtL5U56++uorPPfcc9i2bRv8/f0RERGB3bt3Izk5Gc7Ozli9ejUuXLiA48eP\nAwBKS0vh4eGBgIAAvPvuu7hy5Qqef/55fPDBBwgNDQXQ8NjYMWPGYPny5XjppZdw5swZvPLKK/jy\nyy8xe/ZsAEBmZia8vLzw3HPP4ZVXXkFaWhqWL1+OP/zhD/jrX/+qVWPTIU892V3UFWcvfY/oCwfh\n/8hMBCrn6rqcASkuLg4AoFQqdVwJ9WW8T6g9eJ9QW3iPUHv01s+wOp9DsWDBAhQWFiI8PBw5OTnw\n8vLC0aNHhTUoVCoV0tLShPbm5uaIjo7GihUr4OvrCysrK6xatUoIEwDg5uaGI0eOICwsDJGRkXB0\ndMTWrVuFMAEAzs7OiIqKwsqVK+Ht7Q17e3ssX74c7777bu9dfDc6cDICAPDd2b2YMHo6zOR9M/gQ\nERERkX7ReaAAgJCQEISEhLT43u7du5vt8/T0xKlTpx56zsmTJwvpvTXjx49vc1Xu/uh23g2Mdhun\n6zKIiIiIaADQ6RwK6hlZBem6LoGIiIiIBggGCj2UzUBBRERERL2EgUIPMVAQERERUW9hoNBDOXdu\nQV2v1nUZRERERDQAMFDoCYlYe359XlGWjiohIiIiooGEgUJPaKC9nEhWflorLYmIiIiIug8Dhb54\nYH3CrAIGCiIiIiLqeQwUeuLB5c6z8tN1UQYRERERDTAMFHqKa1EQERERUW9goNAXDwx5KqssRmlF\nkY6KISIiIqKBgoFCTzw4KRsAbufd0EElRERERDSQMFDosVu513VdAhERERHpOQYKPXYrj4GCiIiI\niHoWA4Ue0GiaD3cCgNu5N1p9j4iIiIioOzBQ6IEH508YGRoDAEori1BSUaiLkoiIiIhogGCg0AdN\neiFEEMHFbpjwmvMoiIiIiKgnMVDoG5EIgxXuwsvbnEdBRERERD2IgUIPPDhLwsXufqC4lctHxxIR\nERFRz2Gg0AdaQ56g1UNxK+86J2YTERERUY9hoNADWpOyRSLYWthDJpUDACqqSlFYlqejyoiIiIhI\n3zFQ6BkRRBCJRBisGC7sS8+5qsOKiIiIiEifMVDogZZGNLk5jBC201VXerEaIiIiIhpIGCj0QtMh\nTw1/uNmPFHal5zBQEBEREVHPYKDQA03nUIjuJQo3h/uBIjM/DbV1Nb1eFxERERHpPwYKPdMYKEyM\nzWBn5QQAUNfX4XYeHx9LRERERN2PgUIftPJU2CFNhj2lcdgTEREREfUABgo9oGlhDgWgPewpPSe1\nFysiIiIiooGCgUIfaJrPoQC0J2anqa5wgTsiIiIi6nYMFPpGdD9QONi4wEgqAwCUVhShqCxfV1UR\nERERkZ7qE4Fi+/btGDp0KGQyGZRKJc6cOfPQ9klJSQgICIBcLoeLiws++uijZm1Onz4NpVIJmUwG\nd3d37NixQ+v9PXv2QCwWQyKRQCwWC9s1Nf3vaUit9TuIxRK4Ke6vR3EjO6V3CiIiIiKiAUPngeLA\ngQMIDQ3Fe++9h4SEBEycOBEzZ85EZmZmi+3LysoQGBgIBwcHxMfHY/Pmzfjb3/6GjRs3Cm3S09Mx\na9YsTJo0CQkJCXj77bfx6quv4t///rfWuUxMTKBSqYSvnJwcSKXSHr3enqDRGvKkbajTaGH7RlZS\nL1VERERERAOFga4L2LhxI5YtW4Zly5YBALZs2YJjx44hIiIC4eHhzdrv27cPVVVV2LNnD6RSKUaN\nGoWUlBRs2LABYWFhAICIiAg4OTlh06ZNAICRI0fi559/xvr16zFnzhzhXCKRCIMGDeqFq+xpTSdl\na0cKdydPYft65uXeKoiIiIiIBgid9lDU1tYiPj4egYGBWvuDgoJw7ty5Fo+JjY3F5MmTtXoSgoOD\nkZ2djYyMDKFNUFCQ1nHBwcGIi4uDWq0W9lVVVcHNzQ0uLi747W9/i4SEhO66NJ0RPdBH4WY/AgYS\nQwBAXnE2SioKdVEWEREREekpnQaKgoICqNVqKBQKrf0KhQIqlarFY1QqVYvtNRqNcExrberq6lBQ\nUACgodfi008/xbfffosvv/wSxsbG8Pf3x40b/W8BOE2rsygAQwMpXO2bzKPISu6NkoiIiIhogND5\nkCddmTBhAiZMmCC89vPzg7e3N7Zu3SoMlWpJXFxcb5TXIXfrqoRttVrdrEYTkZWwfe6XH1Bfatxr\ntQ1EffEeob6H9wm1B+8TagvvEXqY4cOH98rn6LSHwtbWFhKJBLm5uVr7c3NzYW9v3+Ix9vb2LbYX\niUTCMa21MTAwgK2tbYvnFYvF8PHxwbVr1zp7Obqjta7dg9OyAYWFq7CdW5LRGxURERER0QCh0x4K\nQ0ND+Pj4IDo6GvPmzRP2R0dHY/78+S0e4+fnh7fffhs1NTXCPIqoqCg4OjrC1dVVaHP48GGt46Ki\noqBUKiGRSFqtJzExEd7e3g+tWalUtuvaelNFdRkOnG/YNjAwaFZjTe0YnEw5AHV9HUqqCjBytDvM\n5JY6qFS/Nf6WqC/eI9R38D6h9uB9Qm3hPULtUVJS0iufo/PHxq5cuRKfffYZdu3ahdTUVLz++uvI\nyclBSEgIAGD16tWYMWOG0H7RokWQy+VYunQpLl++jEOHDmHdunV44403hDYhISHIyspCWFgYUlNT\nsXPnTuzduxerVq0S2nz44YeIiopCWloaEhMTsWzZMiQlJeFPf/pT7118d2ljBWypoREGK9yF19c5\nj4KIiIiIuonO51AsWLAAhYWFCA8PR05ODry8vHD06FE4OzsDaJhgnZaWJrQ3NzdHdHQ0VqxYAV9f\nX1hZWWHVqlUIDQ0V2ri5ueHIkSMICwtDZGQkHB0dsXXrVsyePVtoU1xcjJdeegkqlQoWFhbw9vZG\nTEwMfHx8eu/iu4lWnBA1H/IENDw+Ni0nFUDDehTewyf2fGFEREREpPd0HiiAhh6Fxh6JB+3evbvZ\nPk9PT5w6deqh55w8efJDJypt2LABGzZs6FCdfdXDFrZr5O7shei4rwFwPQoiIiIi6j46H/JE3ayV\nHoqhDh4Qixr+c2ffyUBZZXFvVkVEREREeoqBQi88fA4FABhJZXCzHym8vnIrsScLIiIiIqIBgoFC\nD2jaeGxsIw/XscJ26q3+vyo4EREREekeA4VeaHsOBQB4uN5/JG7qrQStuRdERERERJ3BQKFvWplD\nAQCD7YZBbmQKACitKELOnVu9VRURERER6SkGCj2gacccCgAQiyUYMfgR4TWHPRERERFRVzFQ6IN2\nzqEAAI/B2sOeiIiIiIi6goFCD2i0E8VDeQx+VNi+kXkZtXU1PVQVEREREQ0EDBR6pq0eCmtzO9hZ\nOQEAatU1uJGV3BtlEREREZGeYqDQAx19WtMoVw57IiIiIqLuwUChF5oOeWpjzBMAj8H316NITo/v\niYKIiIiIaIBgoNADmvZPoQAADHceA0OJFACgKryN/OKcnimMiIiIiPQeA4WeaWsOBQBIDY0wssnk\n7KS0Cz1ZEhERERHpMQYKvdDxFa+9ho4XtpNuMlAQERERUecwUOgBTQfnUACA1xClsH0j6zIqq8u7\nuywiIiIiGgAYKPRBB+dQAIC5iRVc7UcAAOo19ZycTURERESdwkChb9rZQwEAXkN8hW3OoyAiIiKi\nzmCg0AOaTsyhAIAxQ+8HiuT0i6hT13ZXSUREREQ0QDBQ6IGmC9u15ylPjRxsXGFtbgcAqK6p5KrZ\nRERERNRhDBT6pv15AiKRCGOaPO0p4fpPPVAQEREREekzBgo905EeCgB4ZNgEYfvX6z9BXa/u7pKI\niIiISI8xUOiBpkOeOmqY4yiYy60AAGVVJbiRdbm7yiIiIiKiAYCBQi90bg4FAIjFEowd7ie8/uXq\n2W6rioiIiIj0HwOFHtDqn+jAY2MbeQ+fJGwn3OCwJyIiIiJqPwYKPdPxOAEMcfSAhYk1AKCiqhTX\nM5O6tygiIiIi0lsMFPqg6RyKTvRQiEVijB0+UXj9y7Uz3VEVEREREQ0ADBR6oLML2zXVdNhT4vVY\nqNV1XT4nEREREek/Bgo9oNVB0alBT4CbwwhYmtoAACqqy3DldmJ3lEZEREREeo6BQt90Lk9ALBJj\n3IjJwuvzKae6px4iIiIi0mt9IlBs374dQ4cOhUwmg1KpxJkzDx/Dn5SUhICAAMjlcri4uOCjjz5q\n1ub06dNQKpWQyWRwd3fHjh07Wj3fF198AbFYjKeeeqrL16IbnX9sbFPjRwUI25du/IyquxVdKYqI\niIiIBgCdB4oDBw4gNDQU7733HhISEjBx4kTMnDkTmZmZLbYvKytDYGAgHBwcEB8fj82bN+Nvf/sb\nNm7cKLRJT0/HrFmzMGnSJCQkJODtt9/Gq6++in//+9/Nznfz5k28+eabmDJlSo9dY0/rysJ2TTna\nusFp0BAAQK26Br9cO9ct5yUiIiIi/aXzQLFx40YsW7YMy5Ytw8iRI7FlyxY4ODggIiKixfb79u1D\nVVUV9uzZg1GjRmHu3Ll46623sGHDBqFNREQEnJycsGnTJowcORLLly/HkiVLsH79eq1z1dXVYdGi\nRfj4448xZMiQHr3OntWkh6ITT3lqavyoqcL2+ZSTXToXEREREem/DgeK7vptOADU1tYiPj4egYGB\nWvuDgoJw7lzLvx2PjY3F5MmTIZVKhX3BwcHIzs5GRkaG0CYoKEjruODgYMTFxUGtvr9o2zvvvIOh\nQ4fi2Wef7a5L6vd8RkyBWNRwW9zMTkF+cY6OKyIiIiKivsygowe4uLhg8eLFWLx4Mby8vLr04QUF\nBVCr1VAoFFr7FQoFTpw40eIxKpUKLi4uzdprNBqoVCq4urpCpVI1CykKhQJ1dXUoKCiAQqFAVFQU\nDh48iMTEjj3NKC4urkPte0NRRa6wXVVV3eUaHS2HIbPoGgDgm5P/wtjBj3fpfANNX7xHqO/hfULt\nwfuE2sJ7hB5m+PDhvfI5He6hGD9+PDZv3oxHH30U3t7e2LhxI1QqVU/U1mMKCgrw/PPPY8+ePTAz\nM9N1OV3WfX1GDYbajRG2b+Zd6tZeKSIiIiLSLx3uoTh06BBKSkpw4MAB7N+/H//v//0/vPnmm5g+\nfTqWLFmC2bNnQyaTtetctra2kEgkyM3N1dqfm5sLe3v7Fo+xt7dvsb1IJBKOaa2NgYEBbG1tcebM\nGahUKkyfPl34Ybm+vh4AIJVKcfny5VYTnVKpbNe19aas/DT8J6FhWy6Xd7nGR+sewYX071F1twLl\nd4thZmcID9ex3VCpfmv8LVFfvEeo7+B9Qu3B+4TawnuE2qOkpKRXPqdTk7ItLCzw4osv4vTp07h5\n8yY++OADZGZmYvHixVAoFHj++edx8mTbE3oNDQ3h4+OD6Ohorf3R0dHw9/dv8Rg/Pz/ExMSgpqZG\n2BcVFQVHR0e4uroKbR48Z1RUFJRKJSQSCXx9fXHp0iUkJCQgMTERiYmJeOqppzBlyhQkJib2uwna\nGq3HxnadoYFUa3L22aTvu+GsRERERKSPuvyUJ1dXV7zzzjs4duwY5s+fj/LycuzZsweBgYFwdXXF\npk2btCZCP2jlypX47LPPsGvXLqSmpuL1119HTk4OQkJCAACrV6/GjBkzhPaLFi2CXC7H0qVLcfny\nZRw6dAjr1q3DG2+8IbQJCQlBVlYWwsLCkJqaip07d2Lv3r1YtWoVgIbf4o8ePVrry9LSEmZmZhg1\nahQMDDrccaNTWiOSuviUp0YTvYKF7Us3fkZJRWG3nJeIiIiI9EuXfnIuKyvDwYMHsW/fPpw+fRoS\niQS/+93vsGTJEkilUvzjH//AypUrkZSUhJ07d7Z4jgULFqCwsBDh4eHIycmBl5cXjh49CmdnZwAN\nk7DT0tKE9ubm5oiOjsaKFSvg6+sLKysrrFq1CqGhoUIbNzc3HDlyBGFhYYiMjISjoyO2bt2K2bNn\nd+Vy+7Dun+PgYOOCYY6jcSM7GfWaesRePoHg8fO7/XOIiIiIqH/rcKBQq9U4evQo9u3bh++++w5V\nVVXw8fHBpk2b8Mwzz8DGxkZo+5vf/AbvvfcetmzZ0mqgABp6FBp7JB60e/fuZvs8PT1x6tSph9Y5\nefLkDj35oKXP6Y+6slL2gyaOCcaN7GQAwE9JUQhUzoVYLOm28xMRERFR/9fhQGFvb4/CwkI4ODjg\n1VdfxZIlSzBq1KhW23t5eaG8vLxLRdLDaT2FqfvyBMa6++HQ6Z2oqC5DYVk+UjJ+gecQTv4iIiIi\novs6PIciODgYx44dw+3bt/HXv/71oWECAJ5++mnhCUrU87qzh8LQQIrHRk8TXp+5dKzbzk1ERERE\n+qHDgeLFF1/E2LFjIWpl8m9BQQF+/PHHLhdG7deT60Q0nZydnBaPvKLsHvssIiIiIup/Ohwopk6d\n2uyRrE2dOHECU6dObfV96lnd2UMBAHZWjhjt5gOg4fG0Pyb+p1vPT0RERET9W4cDRVu/Da+pqYFY\n3OWn0VKHNJ1D0b2BAgCmej8lbMcmn0TlXc6JISIiIqIG7ZqUXVpaiuLiYuH1nTt3cOvWrWbtioqK\n8K9//QtOTk7dVyG1SWsZih44/wiXR+BgMxg5d26hprYaPyUdx3QffX0ELxERERF1RLu6EjZu3Igh\nQ4ZgyJAhEIlECA0NFV43/Ro3bhy+//57vPzyyz1dNzXRk3MoAEAkEiFg7G+F1z8m/hfq+tYXKyQi\nIiKigaNdPRRBQUEwNTWFRqPBm2++iWeeeQbjxo3TaiMSiWBiYgKlUgkfH58eKZbaoQeGPAGA0uNx\nfHvuc1RUlaKoLB+J13/CuBGTeuSziIiIiKj/aFeg8PPzg5+fHwCgoqICc+fOxZgxY3q0MOqI+z0U\n3T0pu5GhgRSTxjyB789/BQA4HncI3sP9W33aFxERERENDB2ePf3+++8zTPQxmp6eRHHPlEd/A0MD\nKQAgM/8mUjIu9tyHEREREVG/0GYPxd69ewEAzz77LEQikfC6Lc8991zXKqMO6Nk5FI3M5JaY6BWE\n0wkNj46NOn8Qo1zHsZeCiIiIaABrM1AsXboUIpEITz/9NKRSKZYuXdrmSUUiEQOFjvTUkKdG08b9\nDmd+PQZ1fR1u5qTgRnYy3J08e/QziYiIiKjvajNQpKWlAQCkUqnWa+o7mj7lqacDhZXZIPiOCkDs\n5eMAgKjz/wf3OQwURERERANVm4HC1dX1oa+pL2i6sF3Pf9oMn7n4OfkkNJp6pN5KQIbqKlztR/T8\nBxMRERFRn8MlrfVA78yguM/OyhHjhvsLr//70796uQIiIiIi6iva7KFYtmxZh08qEomwa9euThVE\nXdPTQ54aPfHYQly8dlbopbiedZlzKYiIiIgGoDYDxcmTJzv8FB8+9ad3aa2U3Uvfe4W1M3w9Hsf5\nlB8AAP89tx+v/T6c/+2JiIiIBpg2A0V6enovlEFd03RSdu+Z+djTiL8SA3V9HW5kJyP1VgJGuXr3\nYgVEREREpGucQ6EHNL09ieIeGwsF/DxnCK//e24/6jX1uimGiIiIiHSCgULf9PKQo6Dx82EoaXik\n8K286/jl6ple/XwiIiIi0q02hzwNGTIEYrEYqampMDQ0xJAhQ9ocJy8SiXDjxo1uK5La0nvrUDzI\n0tQGj499EsfjDwEAvj37OcYMewxSA6NerYOIiIiIdKPNQPH4449DJBJBLBZrvaa+Q3thu94X6DsP\nscknUF5VgqKyfJy6+C2Cxs/XQSVERERE1NvaDBSfffbZQ18TyYxM8JsJz+CrHyIBANFxX2OC5wyY\nm1jpuDIiIiIi6mmcQ6FvdNR75OcVCAebwQCAu7XVOBLLxe6IiIiIBoJOBYq6ujp8+umnmDt3LsaN\nG4dx48Zh7ty5+PTTT1FXV9fdNVIbtIc86SZQSMQSzJ78vPD6p8snkJWfppNaiIiIiKj3dDhQ5OXl\nwdfXF8uXL8epU6cgEokgEolw6tQpLF++HEqlErm5uT1RK7VCg6YL2+mujlGu3hjlOg4AoNHU46sf\ndvAxskRERER6rsOB4tVXX0VycjJ27dqF/Px8xMfHIz4+Hvn5+di5cydSUlLw2muv9USt1A/MmfI8\nJOKGqTlpOamIvXxCxxURERERUU/qcKA4evQoXnvtNTz//POQSCTCfolEgmXLluGVV17B0aNHu7VI\naj9dDXlqZG/tguk+c4TX357Zg7LKEh1WREREREQ9qcOBwsjICK6urq2+P2TIEBgbG3epKOqYpnMo\ndDUpu6mg8b+HjYUCAFB5txzfnPlMtwURERERUY/pcKBYuHAhvvjiC9TW1jZ7r6amBl988QUWLlzY\noXNu374dQ4cOhUwmg1KpxJkzD19tOSkpCQEBAZDL5XBxccFHH33UrM3p06ehVCohk8ng7u6OHTt2\naL1/8OBB+Pr6wsrKCqampvD29sbevXs7VHffoftJ2U1JDYwwP+Al4fX5lB9wLfOSDisiIiIiop7S\n5joU58+f13o9f/58xMTEwNfXFyEhIXB3dwcAXLt2Tfih/fe//327Czhw4ABCQ0MRGRkJf39/bNu2\nDTNnzkRKSgqcnZ2btS8rK0NgYCACAgIQHx+PlJQULF26FKampggLCwMApKenY9asWVi+fDn279+P\nmJgYvPzyy7Czs8OcOQ3DcWxtbfHnP/8ZHh4eMDQ0xHfffYc//vGPsLOzwxNPPNHu+vsCrR6KPmK0\n2ziMHT4RCdfOAQAOnIjAm3/YyBW0iYiIiPSMSNPGT6NisbjZythajym9996D+9RqdbsKmDBhAsaO\nHYvIyEhh34gRIzB//nyEh4c3ax8REYHVq1cjLy8PUqkUABAeHo7IyEjcvn0bAPDWW2/h8OHDuHLl\ninDcCy+8gOTkZJw9e7bVWnx8fPDEE080+9ySkvtzACwsLNp1Xb3pclocdny7FgAw2nUcQmav0XFF\nDUrKC7H28xW4W1MFAJjq/RTmTFmm46p6TlxcHABAqVTquBLqy3ifUHvwPqG28B6h9uitn2Hb7KHY\nvXt3j314bW0t4uPjsWrVKq39QUFBOHfuXIvHxMbGYvLkyUKYAIDg4GCsWbMGGRkZcHV1RWxsLIKC\ngrSOCw4Oxt69e6FWq7Umkzc6ceIErl69iv/5n//phivToT4wh6KRhak15kx+Hl+e2A4AOPXLd3hk\n2GMY5uSp48qIiIiIqLu0GSiWLFnSYx9eUFAAtVoNhUKhtV+hUODEiZYfN6pSqeDi4tKsvUajgUql\ngqurK1QqFQIDA5u1qaurQ0FBgfB5paWlcHJywt27d2FgYIBt27Y1CyIPavyNQF+SWXhN2C4pKelT\nNRpqrOBoORTZxTehgQaf/mc9nhz7Agwl0rYP7qf60vef+i7eJ9QevE+oLbxH6GGGDx/eK5/TqZWy\n9YWZmRkSExMRFxeH8PBwhIWF4YcfftB1WR2mtbBdHyMSieDn/iQMJQ1zJ8qqi3Ax46SOqyIiIiKi\n7tJmD0VLqqurcejQIcTHx6OkpAT19dqrIYtEIuzatavN89ja2kIikTRbWTs3Nxf29vYtHmNvb99i\ne5FIJBzTWhsDAwPY2tpq1Tl06FAAwCOPPILk5GR8/PHHmDp1aqs198WxikY36/FDSsO2paVln6zR\n2EqD/dFbAQBXcuIw/bEn4eE6VsdVdS+OZ6X24H1C7cH7hNrCe4Tao+kcip7U4UBx+/ZtTJ06FTdv\n3oSlpSVKSkpgbW2NoqIi1NfXw9bWFqampu06l6GhIXx8fBAdHY158+YJ+6OjozF//vwWj/Hz88Pb\nb7+Nmpr+Fyo7AAAgAElEQVQaYR5FVFQUHB0dhfUx/Pz8cPjwYa3joqKioFQqW5w/0ai+vh53795t\nV+19VV94bGxLxo+ahoTrP+FyWsNfgPuiNuPNRRthbmKp48qIiIiIqCs6POTpzTffxJ07d3Du3Dlc\nvXoVGo0GBw4cQEVFBcLDwyGXy3H8+PF2n2/lypX47LPPsGvXLqSmpuL1119HTk4OQkJCAACrV6/G\njBkzhPaLFi2CXC7H0qVLcfnyZRw6dAjr1q3DG2+8IbQJCQlBVlYWwsLCkJqaip07d2Lv3r1ak78/\n/vhjnDhxAmlpaUhNTcXf//537Nu3D88++2xHvyU619cWtmuJSCTCM9NXwEzW8ISB0soi7IvejHpN\nfRtHEhEREVFf1uEeiuPHj+NPf/oTJkyYgMLCQgANP9AaGRlh9erVSElJQVhYGL755pt2nW/BggUo\nLCxEeHg4cnJy4OXlhaNHjwprUKhUKqSlpQntzc3NER0djRUrVggL061atQqhoaFCGzc3Nxw5cgRh\nYWGIjIyEo6Mjtm7ditmzZwttysvL8fLLLyMzMxMymQweHh74/PPPsWDBgo5+S/qAvjuHoilzEyss\nDg5FxOEPAACpGb/gh4vfYrrP7DaOJCIiIqK+qsOBory8HMOGDQMAGBndm2hbVia8P2nSJLz55psd\nOmdISIjQI/Gglh5b6+npiVOnTj30nJMnT37okw8+/vhjfPzxxx2qsz/om/0T941y9cZ0nzk4Ef9v\nAMB35z6Hu9NouNqP0HFlRERERNQZHR7y5OTkhOzsbACAiYkJrK2tkZCQILyfkZEBQ0PD7quQ2qQ9\n4qmvRwrgSb8/CAGivl6N3Uf+hvKqUh1XRURERESd0eFAMWXKFHz//ffC63nz5mH9+vVYu3YtPvzw\nQ2zevBnTp0/v1iKpLU2HPPX9QCGRGGDpE29AJpUDAArL8vHZ0fVQ17dvdXUiIiIi6js6HCjCwsIw\nd+5cVFdXAwDWrVsHPz8/rFmzBn/5y18wfvx4bN68udsLpdZpTcruJ2wsFHg2OEx4ffX2r/jPuc91\nWBERERERdUaH51CMGTMGY8aMEV5bWloiOjoaxcXFkEgkMDMz69YCqWP6fv/EfV5DffHEYwtx7OcD\nAIAT8YfhYueOcSMm6bgyIiIiImqvblsp29LSkmFCR7T6J/rBHIqmnnhsITyH3F+U51/RW3E776YO\nKyIiIiKijuhUoCguLsaaNWswbtw4WFtbw9raGuPGjcOaNWtQVFTU3TVSm+5Hir66sF1rxCIxng0O\nxSALBwBATd1d/OPbtSgqK9BxZURERETUHh0OFNevX8cjjzyCtWvXora2FlOnTsXUqVNRW1uLtWvX\n4pFHHsG1a9d6olZqRX+cQ9GU3MgUy3/7jjBJu6SiEP/4Lhx3a6p0XBkRERERtaXDgeKVV15BSUkJ\nTpw4gUuXLuHrr7/G119/jUuXLuH48eMoLS3Fa6+91hO1Unv0rw4KgYONC5bNegtisQQAkJWfhs+O\n/R31fPITERERUZ/W4UARExOD0NBQTJ06tdl706ZNw+uvv44ff/yxW4qjjutvQ56aGjn4USyYen+B\nw8tpcfi/U//s9z0wRERERPqsw4HC0tISVlZWrb5vZWUFS0vLLhVFHdP0B+7+sLDdw0z0CsQMn7nC\n67OXjuG/P/1LhxURERER0cN0OFAsX74cu3btQllZWbP3SkpKsGvXLrzwwgvdUhy1l379Bv9J/8Xw\nGTlFeB114f9w8uI3OqyIiIiIiFrT5joUX331ldZrDw8PiMVijBgxAkuWLIG7uzsA4Nq1a9i7dy8U\nCgVGjhzZM9VSO/TvHgqg4clPiwNfQ9XdCiSnxwMADsfshtzIFBM8uQo7ERERUV/SZqB4+umnIRKJ\nhGE1Tbc/+eSTZu1zc3OxePFiPPPMM91cKrVGa8iTDuvoThKJAZb95k1sP/wX3MxOAQB8cWIbjKUy\njB0+UcfVEREREVGjNgPFDz/80Bt1UBf054XtHkZqaIQXn3oXW7/+M7Ly06DR1OOzY3/HUoChgoiI\niKiPaDNQPP74471RB3WJfs2haEpuZIo//e59bDn4DvKKs1Ffr8ZnR9fjOc1KjBsxSdflEREREQ14\nbQaK1qjValy8eBHp6ekAADc3N/j4+EAs7tTi29RN+vNjY1tjbmKJV+etxdav32sIFZp67D22ARqN\nBj4jJ+u6PCIiIqIBrVM//R84cACDBw/GhAkTsHDhQixcuBATJkyAi4tLs0nc1PO01mnQvzwBALAw\ntcarv18LhZUzADSEiu834kLqaR1XRkRERDSwdThQfPPNN1i0aBEsLS2xZcsWHD9+HMePH8eWLVtg\nZWWFRYsW4dtvv+2JWqkd9LGHopGFiTVenbcW9tYuAACNph6ff78RpxP+o+PKiIiIiAauDg95Cg8P\nh4+PD3788UcYGxsL+6dNm4bly5dj0qRJWLt2LZ566qluLZRaN5BWkm4Y/vQR/vfQGuTcuQUA+Pr0\nTpRVFmOW3x/6/cJ+RERERP1Nh3sokpKSsHjxYq0w0cjIyAjPPvsskpKSuqU46riB8AO1mdwSr/0+\nHG4O99c7ibpwEF+e2A51vVqHlRERERENPB0OFHK5HAUFBa2+n5+fD7lc3qWiqKMGTg9FIxNjM7wy\n50OMdvMR9v10ORq7/rsOd2urdVgZERER0cDS4UAxffp0bN68GTExMc3eO3v2LLZs2YIZM2Z0S3HU\nPtpzsvW/h6KR1NAILzy5GuNHTRX2Jd08j83/9w6KyloPvURERETUfTocKD755BOYmJggICAAvr6+\nWLx4MRYvXgxfX19MmTIFpqamWLduXU/USq3QDMAeikYSiQEWBb6K6T6zhX2Z+Tfx9wOrcCv3ug4r\nIyIiIhoYOhwoXF1dkZiYiNdffx1lZWX4+uuv8fXXX6OsrAxhYWFISEiAq6trT9RK7TEA5lA8SCwS\n43eTluLp6S9DLJYAAEorirD54DtIuHZOx9URERER6bdOLWw3aNAgbNiwARs2bOjueqgzmox5Gnhx\n4r6JXkGwtbDHrv+uQ9XdCtTW1eDTI59ghnIeZvktguRe2CAiIiKi7tOhHorKykpIJBKEh4f3VD3U\nCVoDngZgD0VTI1wewRsLP8EgS0dh3/G4rxHx77+grLJYh5URERER6acOBQq5XI5BgwbBwsKip+qh\nThm4cyhaYmflhJUL12G06zhh39XMS/jkizeQlpOqw8qIiIiI9E+H51AsWLAAX331Ferr63uiHuqi\ngfSUp4cxMTbDi797DzMnPCN8T0rK72DzwXdx6pfvBtRigEREREQ9qcNzKObMmYOTJ09i0qRJeOGF\nFzB06FDIZLJm7caPH98tBVLbmv5wPMBHPGkRi8SY+dhCuCqGY+/3G1FZXYb6ejUO/bgLqRm/YFHg\nazA3sdR1mURERET9WqfWoUhOTkZsbCz++Mc/Ytq0afDz8xO+JkyYAD8/vw6dc/v27UIwUSqVOHPm\nzEPbJyUlISAgAHK5HC4uLvjoo4+atTl9+jSUSiVkMhnc3d2xY8cOrfd37tyJKVOmwNraGlZWVpg2\nbRrOnj3bobr7JiaKB412G4dVz6yHi90wYV9yxkX8df/ruJwWp8PKiIiIiPq/DvdQ7N69u1sLOHDg\nAEJDQxEZGQl/f39s27YNM2fOREpKCpydnZu1LysrQ2BgIAICAhAfH4+UlBQsXboUpqamCAsLAwCk\np6dj1qxZWL58Ofbv34+YmBi8/PLLsLOzw5w5cwA0BI6nn34a/v7+kMvl2LBhA4KDg5GYmIhhw4Y1\n+9y+jMN32mZjrkDo/L/ivz/tw8mL3wAAyqtKsOPbtZjy6G/w1KQlkBoY6bhKIiIiov6n3YGiuroa\n33zzDXJycmBjY4Mnn3wSDg4OXS5g48aNWLZsGZYtWwYA2LJlC44dO4aIiIgWnya1b98+VFVVYc+e\nPZBKpRg1ahRSUlKwYcMGIVBERETAyckJmzZtAgCMHDkSP//8M9avXy8Eis8//1zrvBERETh8+DCO\nHTuGFStWdPm6dIVzKFpnaGCI2ZOfh8dgb+yL3ozSiiIAwI+JR5B6KxGLZryKoY4eOq6SiIiIqH9p\n15Cn7OxseHl5YdGiRXjnnXfw0ksvYcSIETh16lSXPry2thbx8fEIDAzU2h8UFIRz51pekCw2NhaT\nJ0+GVCoV9gUHByM7OxsZGRlCm6CgIK3jgoODERcXB7Va3eJ57969i+rqalhZWXXlknRCa6VsTqJo\nk4frWLz9h83wGnp/nk9eURY2/99qHDq9CzW1d3VYHREREVH/0q4eivfeew/p6ekICwvDtGnTcP36\ndXz00Ud47bXX8Ouvv3b6wwsKCqBWq6FQKLT2KxQKnDhxosVjVCoVXFxcmrXXaDRQqVRwdXWFSqVq\nFlIUCgXq6upQUFDQ7PMar9HMzAxPPfXUQ2uOi+t7Y+4zcjKE7fz8/D5ZY1/kbR8IU5Et4tKOo66+\nBhpocCrhO8SnnoGf+5Owt+jciu/8/lN78D6h9uB9Qm3hPUIPM3z48F75nHYFiujoaCxZsgTr168X\n9ikUCixatAiZmZktznXoTzZv3ox//vOfOHHiBExNTXVdDvUSkUiEEfbj4Gg1DD9d/y9yim8CAMqq\nixCV9Dnc7cZinNs0GBvKdVwpERERUd/VrkChUqng7++vtW/SpEnQaDS4detWpwOFra0tJBIJcnNz\ntfbn5ubC3t6+xWPs7e1bbC8SiYRjWmtjYGAAW1tbrf2bNm3C+++/j2PHjsHHx6fNmpVKZZtteltl\nYh7ON/wsDDs7uz5ZY1/3+MRp+Dn5JP794y5U1VQCAK7nJSC79AZ+O3Ex/LwCIRY9fIRg42+J+P2n\nh+F9Qu3B+4TawnuE2qOkpKRXPqddcyjUajWMjY219jW+rq6u7vSHGxoawsfHB9HR0Vr7o6OjmwWY\nRn5+foiJiUFNTY2wLyoqCo6OjnB1dRXaPHjOqKgoKJVKSCQSYd+GDRvw/vvv48iRIx1+1G3f0mQd\nCk7K7hSRSIQJntPxzrP/izFN5lZUVpfhwMkIbDzwFm7lXtdhhURERER9U7uf8nTz5k2cP39eeN2Y\neFJTU1scJtTehe1WrlyJ5557Dr6+vvD390dERARycnIQEhICAFi9ejUuXLiA48ePAwAWLVqEDz/8\nEEuXLsW7776LK1euYN26dfjggw+Ec4aEhGDbtm0ICwvDSy+9hDNnzmDv3r348ssvhTZ/+9vf8N57\n72H//v1wd3cXejRkMhnMzc3b+23pE7iwXfexMLXGC799B0k3L+Dr0ztxp7ThvsjIvYa/f7kKE8cE\n4zcTnoaZnAviEREREQEdCBTvv/8+3n///Wb7X331Va3XGo0GIpGo1acpPWjBggUoLCxEeHg4cnJy\n4OXlhaNHjwrDqFQqFdLS0oT25ubmiI6OxooVK+Dr6wsrKyusWrUKoaGhQhs3NzccOXIEYWFhiIyM\nhKOjI7Zu3YrZs2cLbbZv3466ujosXLhQq54lS5bg008/bVftpL+8hvpixOBHcPzCIUTHfw21ug4a\naHD20jHEXTmNQOU8BHj/lmtXEBER0YDXrkDR3YvZPSgkJETokWjPZ3t6erb5yNrJkyc/9MkHTUOK\nfmEXRXeRGhjhN37PwHdUAA6e+idSMi4CAO7WVOE/5/bh7K/H8KT/YviMnNLm/AoiIiIifdWuQLFk\nyZKeroO6QHvIEwNFdxtk6YCQ3/0Zl9Pi8M2ZPcgtygQAFJUX4PPvN+HUL9/htxOfFXrniIiIiAaS\ndg95or5La2E76hEikQheQ30xym0cfkqKxpHYL1Be1TCP6HbeDWw//BfYmbtg7OAAAHziBhEREQ0c\nHKehD5gneo1ELMGkR57An5dEIMj39zCU3F+xPa/0NqKSPsf/HlqDtJxUHVZJRERE1HvYQ6Fn+NjY\n3iEzkuPJiYvhP+YJRF04iJ8uR6O+vuFBBFdv/4qrt3/FKNdxCPKdh2FOnjquloiIiKjnMFDoAa0h\nTxzD36uszGyxcFoIZvjMwb4j23Ez71fhv0dKxkWkZFzEUIdRmKGcC88hSs6xICIiIr3DQKEXmi5s\nR7pgY6GA//DfYoyzPzIrLiP+yo9CsLiZk4J/fBcORxtXzFDOhfeISZCIJW2ckYiIiKh/4BwKPaDh\nHIo+w1xmjeeeCMO7z/0vJnjOgER8P7Nn38nA3u834qM9f8LJi9+g8m65DislIiIi6h4MFHqGQ2r6\nBjsrJyya8QrWLI3EVO+nIDU0Ft4rLM3D4ZjdWLNrOf7vh38gtyhLh5USERERdQ2HPOmFpl0UDBR9\niZWZLeZMWYYg39/jx8QjOJ34X1RWlwEAamqrEfPrEcT8egSjXcdhytgn4eE6lovkERERUb/CQKEH\ntBe202Eh1CoTmTlmTnga033mIO7KaZxO+A9y7twS3k/OuIjkjIuwMVfAzysQj42eBgsTax1WTERE\nRNQ+DBR6gFMo+g+poREmegXBzzMQV2//itMJ/8HltDhhAved0lz859w+HPnpX/AaOh4TvYLgMfhR\niDmJm4iIiPooBgq9wy6K/kAkEmHk4EcxcvCjyC/OQcyvR3E++aQwUbteU49fb8Ti1xuxsDYbhAme\nM+A7KgA25godV05ERESkjYFCHzQd8sRA0e8MsnTA3CnL8OTEPyDxeizOJUXhRtZl4f3Csnwcif0C\nR2K/wDDH0fAdFYCxwydCbmSqw6qJiIiIGjBQ6AEubKcfpAZG8PV4HL4ejyO3MBPnkqJwPuUHVNyb\nxA0AN7KTcSM7GQdP/RNeQ33h6xGA0a7jIJHwf2UiIiLSDf4Uog80XNhO3yisnTFnyjI8OfFZ/Hoj\nFhdSfkDqrQTUa+oBAHXqWiRcO4eEa+dgIjOHt/tEjB0+EcOcPLloHhEREfUqBgqiPszQwBA+IyfD\nZ+RklFYUI/7qj7iQegqZeTeFNhVVpThz6RjOXDoGU5kFHh02AWOHT4S7sxfDBREREfU4Bgo9oLUK\nBYc86S1zE0tM9X4KU72fQs6d27iQegpxqadQXH5HaFNeVYKzSd/jbNL3MJVZ4JFhj2Gs+0QMd/bi\nsCgiIiLqEfwJQy9wYbuBxsHGBU/5P4snJ/4BN7NTkHDtLBKu/4TSiiKhTXlVCc4lReFcUhRkUjlG\nufnAa4gSo918IDfmhG4iIiLqHgwUeoAL2w1cYpEY7k6ecHfyxNzHlyMtOwW/XDuHhOvntMJFVU0l\nLl6NwcWrMRCLxBjqNBpjhoyH11BfDLJ00OEVEBERUX/HQEGkJ8QiMYY5eWKYkyfmPv5HpGWnIuH6\nOVy68TMKy/KFdvWaelzPTML1zCT8O+ZTKKydMdp1HDxcvTHMaTSkBkY6vAoiIiLqbxgo9ID2Stns\noqDGcDEaw5xGY+6UPyK7IANJaReQdPM8MnKvabXNLcxEbmEmfvjlWxhKpBjm7AmPwWMxytUb9tYu\nnJdDRERED8VAoQ+4sB09hEgkgtMgNzgNckPw+PkoqSjE5bR4JN08jyu3ElGrrhHa1qprkJrxC1Iz\nfsHhmN2wNLWBx+Cx8HD1xnDnMTCTW+jwSoiIiKgvYqDQA9oL2+muDuofLEysMdErEBO9AlFTexfX\nMi8h9VYCUjJ+QV5Rllbb4vI7iE0+gdjkEwAAB5vBcHfygruzF9ydPBkwiIiIiIGCaCCTGhrBc4gS\nnkOUAIDC0jykZPyC1FsJuHorEVU1lVrtc+7cQs6dW4j59QiAhoAx3NlLCBmmMvNevwYiIiLSLQYK\nfaDVQcEuCuo8a3M7+I8Jhv+YYKjr1chQXUPqrV9w5VYibuVeh7q+Tqt9Y8D4MfF+wBji4CF8DbJ0\n4BwMIiIiPcdAoQe0hzzxhzfqHhKxBEMdPTDU0QO/mfAM7tZWIz3nCq7de0JUeu5V1NertY5pDBjn\nkqIAACYycyFcDHUYCReFO58iRUREpGcYKPRC00nZRD3DyNAYIwc/ipGDHwUA3K2tRlp2Kq5nJeFa\nZhIycq81CxgVVaVIunkeSTfPAwAkYgM4DxqCIQ4ecLUfjsGK4bC1sGcvBhERUT/GQEFEnWJkaAwP\n17HwcB0LoCFgZKiuIS0nFWk5qUjPuYLKu+Vax6jr65CRe03r0bUyIxMMtnOHi8Idrgp3uNi5w8rM\nliGDiIion+gTgWL79u1Yv349cnJy4OnpiU2bNmHSpEmttk9KSsIrr7yC8+fPw8bGBi+++CL+/Oc/\na7U5ffo03njjDVy+fBlOTk5YtWoVXnrpJeH95ORkvP/++7h48SLS0tLwl7/8BWvWrOmxa+xJGq0R\nT/whjHTDyNAYI1zGYITLGAANC+jlFWUhLTv1Xsi4gtyizGbHVd2twJXbibhyO1HYZyazwGDFcLgo\nhsFVMRxOg4bAwsSa9zcREVEfpPNAceDAAYSGhiIyMhL+/v7Ytm0bZs6ciZSUFDg7OzdrX1ZWhsDA\nQAQEBCA+Ph4pKSlYunQpTE1NERYWBgBIT0/HrFmzsHz5cuzfvx8xMTF4+eWXYWdnhzlz5gAAKisr\nMWTIEMybNw/vvfder15z92u6tB1/4KK+QSwSw97aBfbWLvDzCgTQMAQqXXUVaTlXcCvvOm7lXkdl\ndVmzY8uqSnA5PQ6X0+OEfSYyczjbDrm3psYQONm6QWHlDIlE53+NERERDWg6/5d448aNWLZsGZYt\nWwYA2LJlC44dO4aIiAiEh4c3a79v3z5UVVVhz549kEqlGDVqFFJSUrBhwwYhUERERMDJyQmbNm0C\nAIwcORI///wz1q9fLwQKpVIJpbLhUZktfU5/otFwDgX1DyYyc63H1Go0GhSW5iEj9xpu511HRu51\n3M67gbs1Vc2OragqbdaTIZEYwMF6sBAwnAa5wcl2COTGpr12TURERAOdTgNFbW0t4uPjsWrVKq39\nQUFBOHfuXIvHxMbGYvLkyZBKpcK+4OBgrFmzBhkZGXB1dUVsbCyCgoK0jgsODsbevXuhVqshkUi6\n/2KIqMNEIhFsLBSwsVBg3IiGYY71mnrkF2XfCxk3cDv3BrLupLcYMtTqOmTm30Rm/k2t/eYmVnCw\nHgx7Gxc42AyGvfVg2Ns4Q27EoEFERNTddBooCgoKoFaroVAotPYrFAqcOHGixWNUKhVcXFyatddo\nNFCpVHB1dYVKpUJgYGCzNnV1dSgoKGj2ef1d0wFPfGws9XdikRgKa2corJ0xftRUAA0h405JLrIL\n0pGZn4asgnRk5aehqCy/xXOUVhShtKJIqzcDaFgl3N7G5V7YGAx7axc42LhAZmTS49dFRESkr3Q+\n5Km/iYuLa7tRL8vJyRG2s7KyECfqezUOJH3xHtEfhlAYjoDCYQTGOQB366pQVJGLooo8FFbkoqgi\nF8WV+ajXqFs8uqSiECUVhbhySztoyAxNYS6zhrnMBuYyG1jIbGAus4apkSXE4p7p0eR9Qu3B+4Ta\nwnuEHmb48OG98jk6DRS2traQSCTIzc3V2p+bmwt7e/sWj7G3t2+xvUgkEo5prY2BgQFsbW278Qr6\niqZzKNhDQQOHkYEM9hZusLdwE/bVa+pRXl2E4soCFFfmo7gyHyWV+SiputNq0KiqLUdVbTlyS29p\n7ReJxDAzshSChrnM+l7YsIGxoQmfOkVERAQdBwpDQ0P4+PggOjoa8+bNE/ZHR0dj/vz5LR7j5+eH\nt99+GzU1NcI8iqioKDg6OsLV1VVoc/jwYa3joqKioFQquzx/onEid1+SVX0Zl7Matp2cnftkjQNB\n42+J+P3vm9T1ahSUqKC6t5q3qvA2VHduI7c4C2p1XYvHaDT1KK0uRGl1IVB0Tes9Y6kcthb2sLFQ\nwNbCXviysVDAymwQJK30bPA+ofbgfUJt4T1C7VFSUtIrn6PzIU8rV67Ec889B19fX/j7+yMiIgI5\nOTkICQkBAKxevRoXLlzA8ePHAQCLFi3Chx9+iKVLl+Ldd9/FlStXsG7dOnzwwQfCOUNCQrBt2zaE\nhYXhpZdewpkzZ7B37158+eWXQpva2lokJydDo9GguroaKpUKiYmJMDU1xbBhw3r3m9BlfMoTUVsk\nYgkUVk5QWDnhUXc/Yb+6Xo3C0jzkFWUhrzgbeUXZyC/KQm5xNkrK77R6vuqayhYnhAOAWCyBtdmg\ne2HDAbb3QoeNhQK1dXdhaGDUI9dIRESkCzoPFAsWLEBhYSHCw8ORk5MDLy8vHD16VFiDQqVSIS0t\nTWhvbm6O6OhorFixAr6+vrCyssKqVasQGhoqtHFzc8ORI0cQFhaGyMhIODo6YuvWrZg9e7bQJjs7\nG97e3sKQhR07dmDHjh14/PHHcfLkyV66+u7Bhe2IOk8ilmCQpQMGWTrA84H37tZUIa84B/nF2cgt\nykJ+UTbyirKQW5zV4lOnGtXf6w0pKFHhChKbvS81kOHkVQdYmw+CldkgWJvZ3d82t4OJsRn/XyYi\non5DpGm6iAG1qGl3kYWFhQ4radnhmN04efEbAMDvJi3FdJ/ZbRxBPYHdzwOHRqNBWWUJ7pSq7gWH\nXNwpUQkhorSiqEvnlxoaw9psEKzNBsHK3K5h23wQrO4FD3N5z00Wp76Bf59QW3iPUHv01s+wOu+h\nICLqb0QiEcxNLGFuYokhDh7N3q+pvYs7pblCwLhzL3QUlKhQUKxqdXL4/eOrG+Z4FN5u8X2xSAxz\nEytYmtrC0tQGFqbWwrblvW1zE2sYGhh2y/USERE9DAOFHtBaKZujJIh0TmpoBAebwXCwGdzsvQsX\nLqCqphwuQx1QWJqPwrJ8FJXmNfxZlo/C0jzcra1+6PnrNfUoLr+D4ofM8QAAU5nFvZBhA4t7fzZ+\nmZtYw8LECjIjPq2KiIi6hoFCD2iPWeMPBkR9mUgkgtzIDEMcPFrs3dBoNKi8W47C0nwUleU1Cx2F\nZfmoqCpt12eVV5WgvKqkxYnjjQwkhjCXW8LcxLqh10VuBXOTe19Nts3klq0+uYqIiAY2Bgp9oOE6\nFNgYB+gAAB3PSURBVET6QiQSwcTYDCbGZnCxG9pim9q6GpRUFDb0UpQVaG0X39surSiCRlPf5ufV\nqWuFoPLQuiCCicz8XtCwvBc0rIVtU5kFzOQWMJVZwMTYlHM8iIgGEAYKIqJ+xtBAKqx70Rp1vRpl\nlcUPhI4CFJcXorisAKUVRSitLGpzeFUjDTRCj0d2G21FIjFMjM2EgNH4p6nMHGZyyyb7zGEqt4BM\nymFXRET9GQOFHtA0HfTEf5OJCA2Pw22cLwH7Ea22u1tThdLKYpRWFKKkoghllcUNf1YUoaSyCKX3\ntsurSrX/rnkIjaZeCB/tq9UApnILmMnuhwxTY3OYyMxgIvx570tmDhNjMxgaSNt1biIi6nkMFHqG\nQ56IqCOMpDIMksowyNLhoe3U6jqUV5WipKLwXu9GQwgprSwWAkdZVQnKK0tQebe8QzWo6+tQUn7n\noQsJPkhqYKQVMExkZpAbmwlBRG58P4SYyswhNzaDsVTGnhAioh7AQKEHtJ/yxH8siaj7SSQGsDC1\nhoWpdZtt1eo6lFeXoryyBGWVJULQKK+6v930z5p2DrtqqqbuLmrK76KovKD91yA2gNzYtOHLyBQy\nIxPIjUwhNzaBzMhUe9vYFHKj+9tSAyP+/UpE1AoGCiIi6lYSiQEsTKxhYdJ2+AAa1u0or2oIH41/\nVlSXoqKqDBXVTb6qSoXt+vqHr+XREnV9Hcoqi1FWWdzhYyVig4aAcS+MNG63FkqMpSaQGckhk8ph\nLJVDIuE/t0Skv/g3nF7gYudE1H9JDY1gbWgHa3O7drXXaDSorqlqFjoqq8tQfi90VFaX3X/v3r6a\nurudrlFdX4eyez0snWFoIIVMagLjxpBh1BA0GrZNmu3LLs6EVGL8/9u786iozvsN4M+9szCggAur\ngAbcooJapTEuP0XPcYtGY6rWpHUXtwgJaQyuDdQYjbGN2qBpEjUYc1watdYlJ+4LblUbNeKCxn1B\niwdEBGZ9f38MXOc6AyIBxuDzOWfOzH3vd+597/CK83A33M2+CUPRPJ1Gz70kRPRMYqCoBoTjOdn8\nz4aIqjlJkux//ffwKvVKV48zW0xFYSMPBcY85BsfIr8wD/nGPBQUPkS+8dHrAqN62mw1/aI+my0m\nmC0m5OZnP9X7tp5aprzWyFrXgUTvBQ+9Jww6T/uz3hMeqtcGeOiKXhfN413UiagiMVBUA+orrzBQ\nEBG5otPqH1356imZLSZ7wCh8aA8jxUHEMZQYHyLf+BAFhXkoMOWj0PgQhaYCFJjyy3RPkCex2iz2\nvS1lvLFhaTSytiiEGOwhwzGQOIQRvc6gBBTnoOIJD50H9DoDtBod/6BF9BxjoCAiInoCnVYPX23Z\nzwtxJISAyVyIAlM+Coz5KDQVBQ3jQxSW0HY3KxMmixGyFkWh5CGsVkuFbY/VZkF+0aFhFUGWZOh1\nBuh1HvDQGqDXG+zPRYHDQ2coevaAXmuAh94AvdbFPIfp4mXx/BOiZx//lVYHqjtlExHRs0SSJGUv\nQFn3jhw7dgwAEB0drbSZLSYlfNifHR8FMJoKUGguhNFUAKO5wKHN/uz42lYBe0wc2YRN6UtF02i0\nrsNJUSDR6zyKwokHdFo9dFr1tF7rAZ320evH2/U6D2h4Z3eiX4SBohpQnZLNXc5ERNWS/cuyHj41\nav2i5QghYLaaYDQVwlgUMApLCiHmAhSaCl0HE3MhTOZCGM2F5brqVllZrRbkW/Oe+v4mT0OWNfaw\nURw8dEVhQwkoeuh0jtPOAUYdWOy1Wq0eOo0OOq0eWo0eOq0OGlnLw8Oo2mGgqA5Ueyj4S4qIiEom\nSZLy5dkbvhWyTIvVDJPZqAoZJoux6LURJnNB0XNhUV2B/dlin1beUzzfUjxtrJDzT57EZrNW2h6W\nx0mSDJ1GZw8bDoFDp9FDq9VBp7G3O762BxJ13e3M29DIWsgZhU7zdVpdUYBRr0PmnhiqJAwURERE\n9ItoNTpoNTp4GWpW6HKFEEVhxR5MlNBhKXQIMEaYLUaYLCaYLEaYzUaYrSZ7ULEYYXZsL3pdPG0q\nqquK0PJom2xKHyrCgQubylwry5qicOGhCjX2n59W+Tk+emiVoOJ6vuP7SqjROixH82hdDDfVCwNF\nNeB4lSfuRiUioupCkiTlr+w1PCtnHUIIWG2WR2GjKKCop02PwsnjAcZisk87BhhzIcxWMywWE8xW\ns73GaoLFYobVVnEn1z8tm80Ko80KYznuTl/RZEl2HUzKEXA0shYajRZaWQtNUZ1G1irPGk3x68fm\nFS1DNS0z7JQHA0W1wBvbERERlYckScoXU0+PGpW+PqvNCktR2DBZTLAUBw6LCRarCWaLGWZr8fSj\nefZgYoSlaP6t2zdhtVngW8tHFVgcX5usJlWoqco9MU9iq+A9NRVJkuSicOIYVLSPtenU8xzDyxOD\njVYdhB4PNsXrkTWqZcqy5tF8WVNU82yck8NAUQ3wxnZERES/DvYviRp46Az4JfHF1ZXASiOEgM1m\nVYKJ2WIuCjD24GJ/WJymH712bCuhxuKwHKXG9TKepXDzOCFs9lD3C29oWVXkojHlGES0sj1sxA+Y\nWyV9YKAgIiIiquYkSVL+um7QV9LxY0/BcU+N69DxKIw8Cjmuw4nVVvRstcBisz9bbRaHNrN6nsPr\nR21mWIre8yyHHVdsNqs9LMJ9AYiBohrgnbKJiIjo18RxT82zxmazPgocReHkUfgwK23FgUV5bbMH\nHlfBxrndOQQ9WqcZVptVWa/VZoVVaXv02mKzVOolm58GA0V1wBvbEREREVUIWdZAL2sArYe7u/JE\nxYeyWWwWdQApel1VGCiqAd7YjoiIiOj543gomyv379+vkn7IVbIWIiIiIiKqlhgoqgPeKZuIiIiI\n3ISBohpQ39jOjR0hIiIioucOA0W1w0RBRERERFWHgYKIiIiIiMrtmQgUixcvRkREBDw9PREdHY20\ntLRS60+fPo2YmBh4eXkhLCwMs2bNcqrZu3cvoqOj4enpiUaNGuEf//iHU826devQokULGAwGREZG\n4l//+leFbVNVEo7nUPCYJyIiIiKqQm4PFGvWrME777yDGTNm4MSJE+jQoQN69+6NGzduuKx/8OAB\nunfvjuDgYBw/fhwLFy7EJ598gk8//VSpuXLlCvr06YNOnTrhxIkTmDJlCuLi4rBhwwal5tChQxgy\nZAiGDh2KkydP4s0338SgQYNw9OjRSt/miqa+sR0RERERUdVxe6D49NNPMWrUKIwaNQpNmzbFokWL\nEBwcjCVLlrisX7lyJQoKCpCamopmzZrh9ddfR2JiIv72t78pNUuWLEFISAgWLFiApk2bYsyYMRg+\nfDjmz5+v1CxcuBDdunXDlClT0LRpU0ybNg0xMTFYsGBBpW9zhXPIE9xDQURERERVya2Bwmw24/jx\n4+jevbuqvUePHjh48KDL9xw+fBj/93//B71er7T17NkTt27dwtWrV5WaHj16qN7Xs2dPHDt2DFar\n/Rblhw4dcllT0nqJiIiIiMiZW++UnZWVBavVisDAQFV7YGAgdu7c6fI9mZmZCAsLc6oXQiAzMxMN\nGjRAZmamU0gJDAyExWJBVlYWAgMDkZmZ6XK9mZmZpfb5H//+sKybV2Vu3L3kMMU9FERERERUddwa\nKH6NhnSJc3cXnqiqbrNOao0bNwbAz59Kx3FCZcFxQk/CMULPErce8uTn5weNRoM7d+6o2u/cuYOg\noCCX7wkKCnJZL0mS8p6SarRaLfz8/EqtKWm9RERERETkzK2BQqfToW3btti+fbuqffv27ejYsaPL\n97Rv3x779++HyWRS2rZt24Z69eqhQYMGSs3jy9y2bRuio6Oh0WhKrNm+fTs6dOjwi7eLiIiIiOh5\nIQnHmxi4wdq1azFs2DCkpKSgY8eOWLJkCZYvX44zZ84gNDQUU6dOxdGjR7Fjxw4AQG5uLl588UXE\nxMRg+vTpOH/+PEaOHInk5GS88847AOyXjY2KisKYMWMwbtw4pKWlYdKkSVi9ejVee+01APaTsrt0\n6YJZs2bhtddew/r165GUlIQDBw4gOjrabZ8HEREREdGvidsDBQB8/vnnmDdvHm7fvo3IyEgsWLBA\n2UMxcuRI7Nu3Dz///LNSn56ejrfeegv/+c9/ULt2bUyYMAEzZsxQLXP//v1ISEhAeno66tWrhylT\npiA2NlZVs379esyYMQOXLl1Cw4YN8dFHH6F///6Vv8FERERERNXEMxEoiIiIiIjo18ntN7b7NVi8\neDEiIiLg6emJ6OhopKWlubtLVEXmzJmDl156Cb6+vggICEC/fv2Qnp7uVJeUlISQkBB4eXmha9eu\nOHPmjGq+yWRCXFwc/P39UbNmTfTv3x83b96sqs2gKjRnzhzIsoz4+HhVO8cIZWZmYsSIEQgICICn\npyciIyOxf/9+VQ3HyfPNZrNh5syZyneOiIgIzJw5EzabTVXHcfJ82b9/P/r374/Q0FDIsowVK1Y4\n1VTEmMjJycHQoUNRq1Yt1KpVC8OGDSv7VcQElWr16tVCp9OJpUuXinPnzom4uDhRs2ZNcf36dXd3\njapAr169RGpqqkhPTxenT58WAwYMEEFBQSI7O1upmTt3rvDx8REbNmwQ6enpYvDgwaJevXoiLy9P\nqRk/frwICQkRO3fuFD/++KOIiYkRrVu3FjabzR2bRZXk0KFDIjw8XLRu3VrExcUp7RwjlJOTIyIi\nIsSIESPEsWPHxJUrV8SuXbvEuXPnlBqOE5o9e7aoW7eu2LJli7h69arYtGmTqFOnjvjwww+VGo6T\n58/WrVvF9OnTxbp160SNGjVEamqqan5FjYlevXqJyMhIceTIEXH48GHRokUL0a9fvzL1kYHiCdq1\nayfGjRunamvcuLGYNm2am3pE7pSXlyc0Go3YvHmz0hYcHCzmzJmjTBcUFAhvb2/xxRdfCCGEuH//\nvtDr9WLVqlVKzfXr14Usy2Lbtm1V13mqVDk5OaJhw4Ziz549IiYmRhUoOEZo6tSpolOnTqXWcJxQ\n3759xYgRI1Rtw4cPF6+++qoyzXHyfKtZs6ZToKiIMXHmzBkhSZI4dOiQUpOWliYkSRIZGRlP7BcP\neSqF2WzG8ePHne663aNHDxw8eNBNvSJ3ys3Nhc1mQ+3atQEAly9fdrozu8FgQOfOnZUxcuzYMVgs\nFlVNaGgomjVrxnFUjYwdOxaDBw9Gly5dVO0cIwQAGzduRLt27TBkyBAEBgbiN7/5DVJSUpT5HCcE\nAJ06dcLu3btx/vx5AMCZM2ewa9cu9OnTBwDHCTmrqDFx+PBheHt74+WXX1ZqOnbsiBo1apRp3PBO\n2aXIysqC1WpFYGCgqj0wMBA7d+50U6/Ind5++220adMG7du3B2A/JlqSJJdj5NatWwDsN0zUaDSo\nW7euU01mZmbVdJwq1ZdffolLly5h1apVTvM4RggALl26hMWLFyMhIQFTp07FiRMnMGnSJEiShIkT\nJ3KcEAAgMTERDx48QPPmzaHRaGC1WjF9+nSMGzcOAH+fkLOKGhOZmZnw9/d3Wn5AQECZxg0DBVEZ\nvfvuuzh48CAOHDgASZLc3R16RmRkZGD69Ok4cOAAZJk7fck1m82Gl156CbNnzwYAtGrVChkZGUhJ\nScHEiRPd3Dt6VqxevRrffPMNVq9ejebNm+PEiROIj49HeHg4Ro4c6e7uEZWI//uVws/PDxqNBnfu\n3FG137lzB0FBQW7qFblDQkIC1qxZg927dyt3ZAeAoKAgCCFKHSNBQUGwWq24d+9eiTX063Xo0CHc\nu3cPzZs3h06ng06nw969e5GSkgK9Xo+6detyjBCCg4PRrFkzVVuzZs1w7do1APxdQnbvv/8+Jk+e\njEGDBqFFixb4wx/+gHfffRdz5swBwHFCzipqTAQFBeF///uf0/Lv3r1bpnHDQFEKnU6Htm3bYvv2\n7ar27du3Kzfeo+rv7bffVsJE48aNVfPCw8MRFBSkGiOFhYXYv3+/Mkbatm0LrVarqrlx4wbOnj3L\ncVQNDBgwAD/99BNOnjypPKKjo/HGG2/g5MmTaNKkCccIoWPHjspx8cXOnz+v/IGCv0sIAPLz8532\ndMqyrFw2luOEHldRY6J9+/bIy8vD4cOHlZqDBw8iPz8fHTp0eHJHynuW+fNizZo1wsPDQ3z11Vfi\n7NmzIj4+Xnh7e4tr1665u2tUBSZOnCh8fHzE7t27RWZmpvJwvBTbxx9/LGrVqiXWr18vfvrpJ/H7\n3/9ehISEqGomTJggwsLCxI4dO8R///tf0bVrV9GmTRtewq+aevwqTxwjdPToUaHX68Xs2bPFxYsX\nxdq1a4Wvr69YsmSJUsNxQiNGjBBhYWFiy5Yt4sqVK2L9+vXC399fTJ48WanhOHn+5OXliRMnTogf\nf/xReHl5iVmzZokTJ04o30Urakz07t1btGzZUhw6dEgcPHhQREVFif79+5epjwwUZbBkyRIRHh4u\nDAaDiI6OFmlpae7uElURSZKELMtOj+TkZFVdcnKyqFevnvD09BQxMTEiPT1dNd9kMon4+Hjh5+cn\natSoIfr37y9u3LhRlZtCVahr166qQCEExwjZryXfqlUr4enpKZo2bSo+++wzpxqOk+dbXl6eSEhI\nEC+88ILw8vISDRs2FDNmzBBGo1FVx3HyfNmzZ4/L7yMjR45UaipiTOTk5IihQ4cKX19f4evrK4YN\nGybu379fpj5KQghRQXtdiIiIiIjoOcNzKIiIiIiIqNwYKIiIiIiIqNwYKIiIiIiIqNwYKIiIiIiI\nqNwYKIiIiIiIqNwYKIiIiIiIqNwYKIiIiIiIqNwYKIiICACwd+9eyLKMtWvXursrZbZw4UI0atQI\nWq0Wbdq0cVs/RowYgfDwcLetn4jInRgoiIiqUGpqKmRZhsFgwI0bN5zm9+3bFxEREW7omZ0kSW5b\n99M6cOAAEhIS8PLLL2P58uX46KOPSqz9/vvvkZycXGl9kSQJssz/Uono+cTffkREbmA2m11+AXb3\nF3ohhFvX/zT27NkDSZLw+eefY+jQoejVq1eJtVu3bsVf/vKXSuvLV199hXPnzlXa8omInmUMFERE\nbtC6dWssX77c5V6K6i4/P79ClnPnzh0AQM2aNZ9Y+7RBqbCw8KnqNRoNdDrdU72HiKi6YKAgIqpi\nkiRh6tSpEEJg9uzZpdZevXoVsixjxYoVTvNkWVb91T0pKQmyLOPcuXP44x//iFq1asHf3x/Tp08H\nANy8eRMDBgyAr68vgoKCMH/+fJd9s1qt+OCDDxASEoIaNWqgV69euHjxolPthQsXMHjwYPj5+cHT\n0xNt2rTBunXrVDXFh3jt3r0b8fHxCAoKgre3d6nbbLPZMHv2bDRu3BgGgwENGjRAYmKi6ku+LMtI\nSUlRXms0GpefEQCMHDkSixcvVmqL669du6a0TZw4EWvXrkXLli3h4eGhnEfy9ddfo3v37ggODobB\nYECTJk0wd+5cp4Di6hyK4uVu3LgRUVFRMBgMiIyMxA8//KCqe/jwId577z1ERETAYDAgICAAXbt2\nRVpaWqmfExHRs0Lr7g4QET2P6tevj1GjRmHZsmWYPn06QkNDf/Eyiw+XeuONN9CsWTN8/PHH2Lp1\nK+bOnYs6depg2bJl6NKlC+bNm4dvv/0WiYmJiI6ORkxMjLIMIQTmzp0Lm82GyZMnIzs7GwsWLEC3\nbt1w6tQp1KpVCwBw9uxZdOzYEcHBwUhMTETNmjWxYcMGDBo0CCtXrsSbb76p6ltcXBzq1KmDGTNm\n4P79+6Vux9ixY7Fs2TIMHDgQf/rTn3Ds2DF88sknSE9Px+bNmwEAK1euRGpqKnbs2IFvv/0WQgh0\n6NDB5fLGjx+PW7duqWoBwN/fX6nZu3cvvvvuO0yaNAlBQUF48cUXAQCLFy9G8+bN0adPHxgMBuzc\nuRPTpk1Dbm6u6pA1SZJcHq528OBBbNq0CRMmTIC3tzcWLVqEgQMH4tq1a6hdu7bSv+J1N2/eHNnZ\n2Thy5AhOnjyJTp06lfpZERE9EwQREVWZr7/+WsiyLI4cOSKuX78uPDw8xPjx45X5ffv2FeHh4cr0\nlStXhCRJIjU11WlZkiSJ5ORkZTopKUlIkiRiY2OVNqvVKsLCwoRGoxFz5sxR2nNycoSXl5cYOnSo\n0rZnzx4hSZIICgoSubm5SvuuXbuEJEli5syZSlv37t1FVFSUMBqNqj716NFDhIWFqbZXkiTRvn17\nYbVan/j5nDp1SkiSJEaPHq1qT0pKErIsiy1btihtkyZNErIsP3GZT6qVJEloNBpx6tQpp3kFBQVO\nbWPHjhXe3t7CZDIpbSNGjFD93IqX6+HhIS5duuS0fSkpKUpb7dq1RVxcXJm2g4joWcRDnoiI3CQ0\nNBSjRo3C8uXLcf369QpZpiRJGD16tDItyzKio6MhhMCoUaOUdl9fXzRt2hSXLl1yWsbw4cNVhyV1\n7doVLVq0UPYOZGdnY+fOnRg0aBAePHiAe/fuKY+ePXvi5s2buHDhgmqZsbGxZboK0pYtWyBJEhIS\nElTtCQkJkGUZW7ZsKdsH8ZQ6dOiAqKgop3aDwQDAfhhWTk4O7t27h86dO+Phw4dlOgm7W7duqkOh\noqKi4OPjo/rcfX19ceTIEdy6dasCtoSIqOoxUBARudG0adMA4InnUjyN+vXrq6Z9fX2h0+kQEBDg\n1J6dne30/kaNGjm1NWnSBFeuXAEAXLx4EUIIJCUlwd/fX/V47733AAB3795V3itJUpkvhXv16lVI\nkoQmTZqo2n18fBAcHKz0oaI1bNjQZXtaWho6d+4MLy8v1KlTB/7+/hg6dCgAPPHQLQAICwtzaqtd\nu7bqcy8+nKt+/fr47W9/i5kzZyIjI6OcW0JEVPV4DgURkRuFhoZi9OjRWLp0qRIuHJV0GVmbzVbi\nMjUajVNbSXsHRDkuE1u87oSEBLzyyisuayIjI1XTnp6eT72equSqf5cvX0b37t3RtGlTLFy4EGFh\nYTAYDDh+/DimTJlS6s+gmKufBaD+3AcOHIjOnTtj48aN2LZtG/7+979j3rx5SE1NxZAhQ8q/UURE\nVYSBgojIzaZNm4alS5fiww8/dJpXfOJuTk6Oqv3q1auV1p/HD1cCgIyMDLzwwgsAoOxt0Gq16Nat\nW4Wuu0GDBhBCICMjAy1atFDaHzx4gNu3b+PVV18t13LLc3+Pf//73zCZTNi8ebPqpPmff/65XH0o\nTUBAAGJjYxEbG4vc3Fy0a9cOH3zwAQMFEf0q8JAnIiI3CwkJwZgxY5CamuoUFLy9veHn54d9+/ap\n2lNSUirtJngrVqxQHc6za9cupKeno2/fvgDsV0fq2rUrvvzyS5fH/WdlZZV73X369IEQAgsWLFC1\nL1iwADabrdyBokaNGgDKdphSseK9C457IoxGIz777LNy9cEVm82G3NxcVZuPjw/Cw8OdQiQR0bOK\neyiIiKqYq8OMpk6diqVLl+L06dPKnoBiY8aMwdy5cxEbG4vo6Gjs27cPFy5cqLS7WgcGBqJjx44Y\nPXo0srOzsXDhQoSEhKhOlF6yZAk6deqEli1bIjY2Fg0bNsTdu3dx5MgRnD17VnUOwNP0MyoqCqNH\nj8ayZcuQk5ODbt264fjx41i+fDl69+5d6t2wS1N8Yvpbb72F3r17Q6vVol+/fqUeitWzZ0/odDr0\n6dMH48aNQ2FhIVauXAmttuL+63zw4AFCQkLwu9/9Dq1atYKPjw/S0tLwww8/IC4ursLWQ0RUmRgo\niIiqmKs9C8V7KRYvXuw0/89//jOysrLw3Xff4Z///CdeeeUVfP/99wgICCjzXoqS6h5vlyQJiYmJ\nOH/+PObPn4+cnBx07twZixYtUg6/AuwnaR87dgzJycn45ptvkJWVBX9/f7Rq1crp0K2n3ZPyxRdf\nICIiAsuWLcOmTZsQGBiI999/H0lJSWXerse9/vrrSEhIwKpVq7B69WoIIXD58mXUr1+/xHtING7c\nGBs3bsS0adOQmJgIPz8/DB8+HF26dEGPHj2e2JeSluvY7uXlhUmTJmH79u3YtGkTzGYzwsPD8de/\n/hXx8fFl2jYiIneTRGX9iYuIiIiIiKo9nkNBRERERETlxkBBRERERETlxkBBRERERETlxkBBRERE\nRETlxkBBRERERETlxkBBRERERETlxkBBRERERETlxkBBRERERETlxkBBRERERETlxkBBRERERETl\n9v9cnNVbKxM1TwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cb5f212438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from train import Train\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "suite = Train(hypos)\n",
    "suite.Update(60)\n",
    "\n",
    "hs, ps = [], []\n",
    "for hypo, prob in sorted(suite.Items()):\n",
    "    hs.append(hypo)\n",
    "    ps.append(prob)\n",
    "    \n",
    "plt.plot(hs, ps)\n",
    "plt.xlabel('Number of trains')\n",
    "plt.ylabel('Probability');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are too many hypotheses to print, so I plotted the results in\n",
    "Figure [fig.train1]. Not surprisingly, all values of $N$ below 60 have\n",
    "been eliminated.\n",
    "\n",
    "The most likely value, if you had to guess, is 60. That might not seem\n",
    "like a very good guess; after all, what are the chances that you just\n",
    "happened to see the train with the highest number? Nevertheless, if you\n",
    "want to maximize the chance of getting the answer exactly right, you\n",
    "should guess 60.\n",
    "\n",
    "But maybe that’s not the right goal. An alternative is to compute the\n",
    "mean of the posterior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "333.41989326371095\n"
     ]
    }
   ],
   "source": [
    "def Mean(suite):\n",
    "    total = 0\n",
    "    for hypo, prob in suite.Items():\n",
    "        total += hypo * prob\n",
    "    return total\n",
    "\n",
    "print(Mean(suite))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or you could use the very similar method provided by `Pmf`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "333.41989326371095\n"
     ]
    }
   ],
   "source": [
    "print(suite.Mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of the posterior is 333, so that might be a good guess if you\n",
    "wanted to minimize error. If you played this guessing game over and\n",
    "over, using the mean of the posterior as your estimate would minimize\n",
    "the mean squared error over the long run (see\n",
    "<http://en.wikipedia.org/wiki/Minimum_mean_square_error>).\n",
    "\n",
    "You can download this example from <http://thinkbayes.com/train.py>. For\n",
    "more information see Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What about that prior?\n",
    "\n",
    "To make any progress on the locomotive problem we had to make\n",
    "assumptions, and some of them were pretty arbitrary. In particular, we\n",
    "chose a uniform prior from 1 to 1000, without much justification for\n",
    "choosing 1000, or for choosing a uniform distribution.\n",
    "\n",
    "It is not crazy to believe that a railroad company might operate 1000\n",
    "locomotives, but a reasonable person might guess more or fewer. So we\n",
    "might wonder whether the posterior distribution is sensitive to these\n",
    "assumptions. With so little data—only one observation—it probably is.\n",
    "\n",
    "Recall that with a uniform prior from 1 to 1000, the mean of the\n",
    "posterior is 333. With an upper bound of 500, we get a posterior mean of\n",
    "207, and with an upper bound of 2000, the posterior mean is 552.\n",
    "\n",
    "So that’s bad. There are two ways to proceed:\n",
    "\n",
    "-   Get more data.\n",
    "\n",
    "-   Get more background information.\n",
    "\n",
    "With more data, posterior distributions based on different priors tend\n",
    "to converge. For example, suppose that in addition to train 60 we also\n",
    "see trains 30 and 90. We can update the distribution like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for data in [60, 30, 90]:\n",
    "    suite.Update(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these data, the means of the posteriors are\n",
    "\n",
    "      ------- -----------\n",
    "      Upper   Posterior\n",
    "      Bound   Mean\n",
    "      500     152\n",
    "      1000    164\n",
    "      2000    171\n",
    "      ------- -----------\n",
    "\n",
    "So the differences are smaller."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An alternative prior\n",
    "\n",
    "![Posterior distribution based on a power law prior, compared to a\n",
    "uniform prior.](figs/train4.pdf)\n",
    "\n",
    "[fig.train4]\n",
    "\n",
    "If more data are not available, another option is to improve the priors\n",
    "by gathering more background information. It is probably not reasonable\n",
    "to assume that a train-operating company with 1000 locomotives is just\n",
    "as likely as a company with only 1.\n",
    "\n",
    "With some effort, we could probably find a list of companies that\n",
    "operate locomotives in the area of observation. Or we could interview an\n",
    "expert in rail shipping to gather information about the typical size of\n",
    "companies.\n",
    "\n",
    "But even without getting into the specifics of railroad economics, we\n",
    "can make some educated guesses. In most fields, there are many small\n",
    "companies, fewer medium-sized companies, and only one or two very large\n",
    "companies. In fact, the distribution of company sizes tends to follow a\n",
    "power law, as Robert Axtell reports in *Science* (see\n",
    "<http://www.sciencemag.org/content/293/5536/1818.full.pdf>).\n",
    "\n",
    "This law suggests that if there are 1000 companies with fewer than 10\n",
    "locomotives, there might be 100 companies with 100 locomotives, 10\n",
    "companies with 1000, and possibly one company with 10,000 locomotives.\n",
    "\n",
    "Mathematically, a power law means that the number of companies with a\n",
    "given size is inversely proportional to size, or\n",
    "\n",
    "$${\\mathrm{PMF}}(x) \\propto \\left( \\frac{1}{x} \\right)^{\\alpha}$$ \n",
    "\n",
    "where\n",
    "${\\mathrm{PMF}}(x)$ is the probability mass function of $x$ and $\\alpha$\n",
    "is a parameter that is often near 1.\n",
    "\n",
    "We can construct a power law prior like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "class Train(Dice):\n",
    "\n",
    "    def __init__(self, hypos, alpha=1.0):\n",
    "        Pmf.__init__(self)\n",
    "        for hypo in hypos:\n",
    "            self.Set(hypo, hypo**(-alpha))\n",
    "        self.Normalize()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here’s the code that constructs the prior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from train3 import Train\n",
    "\n",
    "hypos = range(1, 1001)\n",
    "suite = Train(hypos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the upper bound is arbitrary, but with a power law prior, the\n",
    "posterior is less sensitive to this choice.\n",
    "\n",
    "Figure [fig.train4] shows the new posterior based on the power law,\n",
    "compared to the posterior based on the uniform prior. Using the\n",
    "background information represented in the power law prior, we can all\n",
    "but eliminate values of $N$ greater than 700.\n",
    "\n",
    "If we start with this prior and observe trains 30, 60, and 90, the means\n",
    "of the posteriors are\n",
    "\n",
    "      ------- -----------\n",
    "      Upper   Posterior\n",
    "      Bound   Mean\n",
    "      500     131\n",
    "      1000    133\n",
    "      2000    134\n",
    "      ------- -----------\n",
    "\n",
    "Now the differences are much smaller. In fact, with an arbitrarily large\n",
    "upper bound, the mean converges on 134.\n",
    "\n",
    "So the power law prior is more realistic, because it is based on general\n",
    "information about the size of companies, and it behaves better in\n",
    "practice.\n",
    "\n",
    "You can download the examples in this section from\n",
    "<http://thinkbayes.com/train3.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credible intervals\n",
    "\n",
    "Once you have computed a posterior distribution, it is often useful to\n",
    "summarize the results with a single point estimate or an interval. For\n",
    "point estimates it is common to use the mean, median, or the value with\n",
    "maximum likelihood.\n",
    "\n",
    "For intervals we usually report two values computed so that there is a\n",
    "90% chance that the unknown value falls between them (or any other\n",
    "probability). These values define a **credible interval**.\n",
    "\n",
    "A simple way to compute a credible interval is to add up the\n",
    "probabilities in the posterior distribution and record the values that\n",
    "correspond to probabilities 5% and 95%. In other words, the 5th and 95th\n",
    "percentiles.\n",
    "\n",
    "`thinkbayes` provides a function that computes percentiles:\n",
    "\n",
    "```python\n",
    "def Percentile(pmf, percentage):\n",
    "    p = percentage / 100.0\n",
    "    total = 0\n",
    "    for val, prob in pmf.Items():\n",
    "        total += prob\n",
    "        if total >= p:\n",
    "            return val    \n",
    "```\n",
    "\n",
    "And here’s the code that uses it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(50, 950)\n"
     ]
    }
   ],
   "source": [
    "from thinkbayes import Percentile\n",
    "\n",
    "interval = Percentile(suite, 5), Percentile(suite, 95)\n",
    "print(interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the previous example—the locomotive problem with a power law prior\n",
    "and three trains—the 90% credible interval is $(91, 243)$. The width of\n",
    "this range suggests, correctly, that we are still quite uncertain about\n",
    "how many locomotives there are."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cumulative distribution functions\n",
    "\n",
    "In the previous section we computed percentiles by iterating through the\n",
    "values and probabilities in a Pmf. If we need to compute more than a few\n",
    "percentiles, it is more efficient to use a cumulative distribution\n",
    "function, or Cdf.\n",
    "\n",
    "Cdfs and Pmfs are equivalent in the sense that they contain the same\n",
    "information about the distribution, and you can always convert from one\n",
    "to the other. The advantage of the Cdf is that you can compute\n",
    "percentiles more efficiently.\n",
    "\n",
    "`thinkbayes` provides a `Cdf` class that\n",
    "represents a cumulative distribution function. `Pmf` provides\n",
    "a method that makes the corresponding Cdf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cdf = suite.MakeCdf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And `Cdf` provides a function named `Percentile`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(50, 950)\n"
     ]
    }
   ],
   "source": [
    "interval = cdf.Percentile(5), cdf.Percentile(95)\n",
    "print(interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Converting from a Pmf to a Cdf takes time proportional to the number of\n",
    "values, `len(pmf)`. The Cdf stores the values and\n",
    "probabilities in sorted lists, so looking up a probability to get the\n",
    "corresponding value takes “log time”: that is, time proportional to the\n",
    "logarithm of the number of values. Looking up a value to get the\n",
    "corresponding probability is also logarithmic, so Cdfs are efficient for\n",
    "many calculations.\n",
    "\n",
    "The examples in this section are in <http://thinkbayes.com/train3.py>.\n",
    "For more information see Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The German tank problem\n",
    "\n",
    "During World War II, the Economic Warfare Division of the American\n",
    "Embassy in London used statistical analysis to estimate German\n",
    "production of tanks and other equipment.[^2]\n",
    "\n",
    "The Western Allies had captured log books, inventories, and repair\n",
    "records that included chassis and engine serial numbers for individual\n",
    "tanks.\n",
    "\n",
    "Analysis of these records indicated that serial numbers were allocated\n",
    "by manufacturer and tank type in blocks of 100 numbers, that numbers in\n",
    "each block were used sequentially, and that not all numbers in each\n",
    "block were used. So the problem of estimating German tank production\n",
    "could be reduced, within each block of 100 numbers, to a form of the\n",
    "locomotive problem.\n",
    "\n",
    "Based on this insight, American and British analysts produced estimates\n",
    "substantially lower than estimates from other forms of intelligence. And\n",
    "after the war, records indicated that they were substantially more\n",
    "accurate.\n",
    "\n",
    "They performed similar analyses for tires, trucks, rockets, and other\n",
    "equipment, yielding accurate and actionable economic intelligence.\n",
    "\n",
    "The German tank problem is historically interesting; it is also a nice\n",
    "example of real-world application of statistical estimation. So far many\n",
    "of the examples in this book have been toy problems, but it will not be\n",
    "long before we start solving real problems. I think it is an advantage\n",
    "of Bayesian analysis, especially with the computational approach we are\n",
    "taking, that it provides such a short path from a basic introduction to\n",
    "the research frontier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "Among Bayesians, there are two approaches to choosing prior\n",
    "distributions. Some recommend choosing the prior that best represents\n",
    "background information about the problem; in that case the prior is said\n",
    "to be **informative**. The problem with using an\n",
    "informative prior is that people might use different background\n",
    "information (or interpret it differently). So informative priors often\n",
    "seem subjective.\n",
    "\n",
    "The alternative is a so-called **uninformative prior**,\n",
    "which is intended to be as unrestricted as possible, in order to let the\n",
    "data speak for themselves. In some cases you can identify a unique prior\n",
    "that has some desirable property, like representing minimal prior\n",
    "information about the estimated quantity.\n",
    "\n",
    "Uninformative priors are appealing because they seem more objective. But\n",
    "I am generally in favor of using informative priors. Why? First,\n",
    "Bayesian analysis is always based on modeling decisions. Choosing the\n",
    "prior is one of those decisions, but it is not the only one, and it\n",
    "might not even be the most subjective. So even if an uninformative prior\n",
    "is more objective, the entire analysis is still subjective.\n",
    "\n",
    "Also, for most practical problems, you are likely to be in one of two\n",
    "regimes: either you have a lot of data or not very much. If you have a\n",
    "lot of data, the choice of the prior doesn’t matter very much;\n",
    "informative and uninformative priors yield almost the same results.\n",
    "We’ll see an example like this in the next chapter.\n",
    "\n",
    "But if, as in the locomotive problem, you don’t have much data, using\n",
    "relevant background information (like the power law distribution) makes\n",
    "a big difference.\n",
    "\n",
    "And if, as in the German tank problem, you have to make life-and-death\n",
    "decisions based on your results, you should probably use all of the\n",
    "information at your disposal, rather than maintaining the illusion of\n",
    "objectivity by pretending to know less than you do."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "To write a likelihood function for the locomotive problem, we had to\n",
    "answer this question: “If the railroad has $N$ locomotives, what is the\n",
    "probability that we see number 60?”\n",
    "\n",
    "The answer depends on what sampling process we use when we observe the\n",
    "locomotive. In this chapter, I resolved the ambiguity by specifying that\n",
    "there is only one train-operating company (or only one that we care\n",
    "about).\n",
    "\n",
    "But suppose instead that there are many companies with different numbers\n",
    "of trains. And suppose that you are equally likely to see any train\n",
    "operated by any company. In that case, the likelihood function is\n",
    "different because you are more likely to see a train operated by a large\n",
    "company.\n",
    "\n",
    "As an exercise, implement the likelihood function for this variation of\n",
    "the locomotive problem, and compare the results."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
